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A QUMOND galactic N-body code I: Poisson solver and rotation 
curve fitting 
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ABSTRACT 

Here we present a new particle-mesh galactic N-body code that uses the full multigrid algo- 
rithm for solving the modified Poisson equation of the Quasi Linear formulation of Modified 
Newtonian Dynamics (QUMOND). A novel approach for handling the boundary conditions 
using a refinement strategy is implemented and the accuracy of the code is compared with 
analytical solutions of Kuzmin disks. We then employ the code to compute the predicted rota- 
tion curves for a sample of five spiral galaxies from the THINGS sample. We generated static 
N-body realisations of the galaxies according to their stellar and gaseous surface densities and 
allowed their distances, mass-to-light ratios (M/ L) and both the stellar and gas scale-heights 
to vary in order to estimate the best fit parameters. We found that NGC 3621, NGC 3521 and 
DDO 154 are well fit by MOND using expected values of the distance and M / L. NGC 2403 
required a moderately larger M/L than expected and NGC 2903 required a substantially 
larger value. The surprising result was that the scale-height of the dominant baryonic compo- 
nent was well constrained by the rotation curves: the gas scale-height for DDO 154 and the 
stellar scale-height for the others. In fact, if the suggested stellar scale-height (one-fifth the 
stellar scale-length) was used in the case of NGC 3621 and NGC 3521 it would not be pos- 
sible to produce a good fit to the inner rotation curve. For each of the four stellar dominated 
galaxies, we calculated the vertical velocity dispersions which we found to be, on the whole, 
quite typical compared with observed stellar vertical velocity dispersions of face on spirals. 
We conclude that modelling the gas scale-heights of the gas rich dwarf spiral galaxies will be 
vital in order to make precise conclusions about MOND. 



1 INTRODUCTION 



The H I Nearby Galaxy Survey (THINGS . Iwalter et alj2o63) 
brought an unprecedented level of precision to the measurement of 
the rotation curves of certain nearby spiral galaxies. This, when 
coupled with 3.6 /im images of the stellar com ponent from the 
Spitz er Infrared Nearby Galaxies Survey (SINGS. lKennicutt et al.l 
l2003h . produces a stringent new data set for modeling the galac- 
tic dynamics of these systems. These tighter constraints are of 
paramount importance for testing alternative theories of grav- 
ity, in particular those with no gala ctic dark matte r like Modi- 
fied Newtonian Dynamics (MOND - iMilgroml [19831 but see also 
iFamaev & McGaughll2oTTh . 

MOND is an appealing framework because, at least in galax- 
ies, the gravitational field is fully determined by the matter distribu- 
tion of the visible components. This means that a galaxy comprised 
of disky stellar and gaseous components, produces a rotation curve 
depending only on the properties of those components. It is cru- 
cial in order to keep up with the advancing observations, that we 
produce methods of modelling the galaxies with a similar level of 
sophistication. 



In this article, we develop a galactic N-body code for MOND 
and apply it to fitting the high resolution rotation curves of the 
THINGS survey. We compare predicted and fitted distances and 
make use of free parameters in the form of the scale-heights of the 
two baryonic components. 



2 SOLVING THE MOND EQUATION ON A GRID 

MOND has a different force law from Newtonian 
dynamics, and this traditionally is rather tricky to solve 
jBekenstein & Milgromlll984 IMilgroml 1 198d liirada & Milgroml 
1 1995b . For instance, in the Newtonian analogue, the ordinary Pois- 
son equation must simply be solved using the matter distribution 
in terms of its density, which includes stars (p„), gas (p s ) and cold 
dark matter (pcdm), using 



V $]v = AitG{p, + p g + pcdm) 



(1) 



to give the Newtonian potential, $jv. Howeve r, in a recently pro- 
posed version of MOND, called QUMOND dMilgrordlioTcl. but 
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see also lZhao & Famaevll2010l) it was shown that the MOND po- 
tential, $, can be found exactly from the Newtonian potential (not 
including cold dark matter) as follows 

V 2 $ = V • {v(y)V® N ] , (2) 

where u(y) = 0.5 + 0.5^1 + 4/y and y = V$jv/a D , with 
a being the MOND acceleration constant chosen here to be 
3.6 (kms -1 ) pc~ . v(y) could take on another form, and simi- 
larly a a could take on a different value, but we choose here not to 
focus on this topic for fear of being sidetracked. We note that the 
i/-function we adopt here is the analogue of the /i-function used by 
iFamaev & BinnevT ( l2005l) to fit the terminal velocity curve of the 
Milky Way. 

Very often in the literature, Eq[2]is not solved, but rather an 
approximation is made that is referred to as the algebraic MOND 
relation 

V$ = v{y)V$> N . (3) 

This approximation ignores the curl field, which is negligible in 
most cases of high symmetry, but has the scope to be a significant 
component of the gravitational field in triaxial systems. Here we are 
not advocating that it is imperative to use Eq[2]for fitting rotation 
curves of relaxed spiral galaxies, but rather that when MOND is 
extended to study interacting galaxies, or galaxies with lopsided 
matter distributions. 

Our goal is to determine the gravity at locations inside the 
galaxy we are studying and the procedure to do this is as follows. 
We must set up a three dimensional grid, sliced into numerous cubic 
cells - usually we use 65, 129 or 257 cells per dimension. In each 
cell there are the following. 

• An approximated, but near exact, value of the baryonic density 
(step 1, the source of Eq[Tjusing only gas and stars). 

• A solution to the Newtonian potential (step 2). 

• The source "density" - the right hand side of Eq[2](step 3). 

• The solution of the QUMOND potential (step 4) 

Getting from step 1 to 4 is in the details of the code, which we 
briefly describe here. 

We read in a set of particle positions that determine the density 
of the stars and gas in the galaxy on the grid using the cloud-in-cell 
technique (step 1). The stars and gas are allocated half of the total 
number of particles each and therefore the particles of gas and stars 
are weighted differently and according to the total mass of each 
component. With this density, we solve for the Newtonian potential 
(step 2; EqQJ usi ng the full multigr id algorithm that is well known 
and described in IPress etalJ l ll992h . The full multigrid method is 
extended to three dimensions which implements red-black Gauss- 
Seidel relaxation. Our default is to use six V-cycles, from finest 
to coarsest grid and back to finest with bilinear prolongation and 
restriction, as well as two pre and post smoothing steps before and 
after the coarse grid correction is computed. 

We then use finite differencing to find the source for Eq[2](step 
3) . The finite differencing te chnique used to solve Eq[2]is described 
in I Angus & Diafericl J201 lh . We use the same multigrid technique 
to solve EqQ]and Eq[2] only the source changes (step 4). We make 
a final finite differencing to determine the QUMOND gravity, V$, 
on the grid and then we interpolate to the point at which we wish 
to know the gravity. The interpolation scheme is the same as the 
scheme to assign density to the grid i.e. the cloud-in-cell technique. 

Our final step is to calculate the circular speed associated with 



this QUMOND potential as a function of radius. This is taken as 

vl = R x 0h*|*=o. (4) 

Although this is the first attempt to solve the modified Pois- 
son equat ion of QUMOND on a grid in a no n -cosmological set- 
ting (see iLlinares etail 120081 ; iLlinaresI l201ll : I Angus & Diafericl 
1201 ll for the cosmological analogue), this is not the first time a 
modified P oisson equation linked to MOND has been solved in 
such a way. iBrada & Milgroml dl999h developed a cubic particle- 
mesh code that solved the modified Poisson equation and used it 
to investigate some imp ortant topics jBrada & Milgromll2000b1 lal). 
iTiret & Combes! | |2007|) developed a similar code and parallelised 
it and incorporated hy drodynamics into the ir simulations (see 
ITiret & Combesll2008al fbl). In addition to these. lLondrillo & NipotH 
d2009h produced a spherical grid code that investigat ed variou s is- 
sues in the framework of MOND dNipoti et alj2007al lbl l200l) . 

The key to a galactic N-body code, as opposed to a cosmolog- 
ical one, is the boundary conditions. In a cosmological code, the 
expectation that the Universe is homogeneous and isotropic allows 
us to enforce periodic boundary conditions. At the galaxy scale, 
the boundary conditions must be set precisely at all outer grid cells, 
which is non-trivial in MOND. It is not clear how accurately this 
was achieved in the previous cod es, but here we implemen t a differ- 
ent strategy (see lWu et alj2009l and lTiret & Combesl20"o"7l for other 
examples). 

We define a coarse grid which is many times (2 10 ) larger than 
the galaxy we are studying. For this coarsest grid we set the bound- 
ary condition of the Newtonian and QUMOND potentials to be 
zero. Then we solve for the Newtonian and QUMOND potentials 
everywhere on the grid. From this, we define a refined grid that is 
half the size of the coarse grid and we interpolate through the values 
of the potential on the coarse grid to define the boundary condition 
of the refined grid. We make this refinement up to ten times in order 
to zoom in on our galaxy, but the boundary conditions are correct 
at the sub percent level at all points on the grid by the second re- 
finement. This also enables us to use different grid levels to find 
the potential at the location of a particle, depending on its position 
in the various grid levels. The limitation of this is that it is cen- 
trally refining, not arbitrarily towards regions of high density. Note 
that the number of grid cells is constant for each refinement, so the 
resolution increases for the finer grids - with more cells per kpc. 

There is an initial guess that is required in multigrid methods. 
For clarity, our finest grid is 257x257x257 and we solve for the po- 
tential on this grid at 10 incremental levels of resolution (box size). 
For each box size, in order to more rapidly calculate the potential 
on the finest grid, we smooth the calculated (fed in) density down 
to the coarsest possible level of resolution, which is a 3x3x3 grid. 
The coarsest density field is known in each of the 27 cells and we 
must solve for the potential in each of these cells. Since for the 
largest box, of length 4 Mpc, the density and potential in each of 
the 26 outer cells is zero, the potential in the central cell is sim- 
ply proportional to the density in that central cell. This is true for 
both the Newtonian step and the QUMOND step. Following this, 
we perform the standard multigrid interpolation of the potential to 
finer grids and use the finer density defined on those grids to make 
a more rapid calculation of the potential on the finest grid. 



3 STATIC N-BODY REALISATIONS 

Our aim is to use this particle mesh code to fit rotation curves. 
To do this, we need to generate static N-body realisations of galax- 
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ies depending on their surface density. Initially, to test the resolu- 
tion of the code, we generated a realisation of a Kuzmin disk and 
compared the analytical circular velocity to the circular velocity 
simulated by the code. 



3.1 Kuzmin Disks 

Kuzmin disks are unique in MOND because their isopoten- 
tial contours are spherically symmetric (at least in each distinct 
hemisphere, above and below the disk mid-plane), meaning there 
is no curl field. The surfac e density of a Kuzmin disk is given by 
iBinnev & Tremaind d2008h to be E K (R) = ^(R 2 + a 2 )~ 3/2 . 
Assuming an infinitely thin disk, and integrating this by 2nRdR, 
we get the enclosed mass M K (R) = M(l r =2 ). If 

V R 2 +a 2 

we express the fraction of the total mass as / = Mk{R)/M, 
then we know the way to distribute particles in radius is R — 
a\J (1 — /)~ 2 — 1. From this we randomly sample / and distribute 
the particles according to a random azimuthal angle. 

We did this for a model where we set a = 1.5 kpc and 
M = 10 10 M© and we show in Fig Q] (left hand panel) the frac- 
tional difference between this, the simulated QUMOND rotation 
curve, and the analytically expected rotation curve. For the Kuzmin 
disk, the 3D Newtonian gravity on a star slightly below the disk 
is defined by gN,r = —GM/(R + a 2 ) and this points towards 
the position (R, z) — (0, a), as if all the mass was concentrated 
in a single point at the location (0, a). Similarly, the gravity to- 
wards the centre of the disk is g N , R = —GMR(R 2 + a 2 )~ 3/2 . 
Therefore, the QUMOND gravity towards the centre of the disk is 
9m,r = v(\gN tT \/a )gN,R- Note the subtlety here with respect to 
the subscripts of R and r. 

Fig demonstrates that the accuracy of the code is less than 
1% beyond 0.3 kpc and less than 0.1% beyond 1 kpc for this model. 
This depends greatly on the extent of the galaxy. A Kuzmin disk, 
with scale radius a — 1.5 kpc still has 5% of its mass beyond 
30 kpc and so is rather unrealistically extended. In Fig[T](right hand 
panel) we make the same plots for a — 150 pc and M — 1O 9 M0 
to demonstrate that in that case the agreement is better than 0. 1 % 
down to ~50 pc. In both cases, as for the real galaxies to follow, 
we use the following strategy to evaluate gravity at various radii, R. 
When R < 0.8 kpc we use the finest grid which has a default box- 
size of 4 kpc, when 0.8 < R < 1.6 kpc we use the second finest 
box of size 8 kpc, for 1.6 < R < 3.2 kpc we use the 16 kpc box 
and so on. These values need not be fixed and should be tailored to 
different galaxy sizes. 



3.2 Convergence 

Although we can accurately simulate the gravitational field 
of a Kuzmin disk to roughly 0.1%, it is equally important to en- 
sure that the solution to Poisson's equation does converge for more 
complex density distributions, which might not immediately follow 
from a smooth Kuzmin disk. As mentioned above, our default pa- 
rameters for using the full multigrid algorithm are that we make a 
standard number of six V-cycles at each level of the grid. Before 
and after we make the coarse grid correction we always make two 
red-black Gauss-Seidel relaxation sweeps. 

For one of our galaxies, NGC 3621, we have plotted in Fig 
[3] the relative error (as a function of radius) between the simulated 
rotation speed for a series of numbers of V-cycles compared to the 
rotation speed found using 99 V-cycles. Each different number of 



V-cycles is given a distinct colour. If a colour is not seen at a spe- 
cific radius, this is because the relative error is less than 10 -8 and 
is most likely identical to the 99 V-cycles simulation. The discon- 
tinuities are caused by moving from one grid resolution to another. 
Clearly, using only one or two V-cycles does not allow for conver- 
gence, but by six cycles the difference is less than 0. 1% at all radii. 
Since the theoretical accuracy of the code is roughly 0.1% as well 
(see Fig|TJ, it makes 6 V-cycles a sensible number of cycles to make 
until the accuracy can be improved. 

3.3 Curl-field and scale-height 

To compare the difference between solving the modified Pois- 
son equation of QUMOND (EqO and the algebraic MOND re- 
lation (Eq[3) we plot in Fig [2] the rotation curve for an exponen- 
tial disk galaxy with scale length a = 1.5 kpc and total mass 
M = 1O 1O M . There is a small but noticeable difference even 
for a perfectly symmetric and thin exponential disk. Also plotted is 
the simulated rotation curve for a disk with scale-height of 1 kpc, 
with a sech 2 distribution. This creates a 10 kms -1 difference out 
to R — 6 kpc and is still significant at R = 10 kpc. For this reason, 
scale-height is a crucial parameter that we will use to enhance the 
fits to the rotation curves of the galaxies. 

3.4 Realistic two-component galaxies 

In order to generate realistic, static N-body realisations for 
galaxies we need the surface densities of both the stellar disk and 
the gaseous disk, like that plotted for NGC 3621 in Figg] In our 
models we always separate stars and gas, giving half of the particles 
to each component and weighting them according to the relative 
masses of thes e two components. We use the well known rejection 
technique from lPress et al.l 11992) to produce an N-body realisation 
of the two components that resembles the observed surface densi- 
ties. Typically we use 256 3 particles and our most refined grid is 
only 4 kpc across, with 129 cells per dimension. 

3.5 Free parameters 

We are now in a position to compare simulated circular ve- 
locity curves with the observed rotation curves. This brings into 
question the various free parameters we employ. We assume that 
the uncertainty in the inclination of the galaxy is contained within 
the measurement errors associated with each data point. The er- 
ror from inclination is addressed in de Blok et al. (2008), where 
the tilted ring fits find little variation in inclination for the galax- 
ies we use. This leaves us with four free parameters: the distance 
to the galaxy, D, the mass-to-light ratio (M / L) of the stellar com- 
ponent of the galaxy and the scale-heights of both the stellar, z*, 
and gas, z g , disks - with sech 2 distributions. We then make an ex- 
haustive search through parameter space to find the quality of fit 
values, X 2 /n = (e? =1 W*0) 2 \ / n> f these param- 

eters by comparing the simulated QUMOND rotation curve with 
the observed one. We initially enforce flat priors on the distance to 
be no more than two standard deviations and the M/L is forced to 
be within the range set by the diet-Salpeter and Kroupa IMFs. We 
then relax these priors if no good fit is found, or if a significantly 
improved fit is found despite not being preferred by current limits 
on the parameters. The scale-heights are free to vary from razo r thin 
up to 1 kpc. It should be emphasise d that lGentile et al.l d201 lh have 
suggested that lde Blok et al .1 d2008l) overestimate the error bars on 
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some rotation curves, a a t s , and as such the reduced \ 2 cannot be 
directly used as a probability indicator in every case. 



4 SAMPLE 

The sample of five galaxies that we've chosen is a subsam - 
ple of the galaxies used in the MOND fits of lGentile et ail j201 lh. 
which is itself a subset of the galaxies studied by Ide Blok et al .1 
d2008l) . Our sample removed all galaxies that had a significant 
spheroid component and also those galaxies that could be aided 
by the external field effect, which we devote a future paper to. Here 
our goal is not to have a representative sample of galaxies to test 
MOND in a statistical sense, but rather to highlight the potential of 
this Poisson solver and to investigate the impact of variable scale- 
heights on rotation curve fits. The properties of t he fiv e galaxies 
themselves can be found detailed in Ide Blok et alj<2008l) . 



5.1 DDO 154 

A lower \ 2 can t> e achieved by increasing the distance and gas 
scale-height. In table[T]we give the parameters used for the three fits 
in the left hand panel of Fig [6] In all figures showing the rotation 
curve, the fits and data points are rescaled to the distance of model 
(a) for that particular galaxy. In the case of DDO 154 this is 1.085 
times the measured distance of 4.3 Mpc and as such is 4.67 Mpc. It 
is not possible to fit the inner curve and outer flat data points beyond 
7 kpc with the same model. The flat ou ter points may be due to a 
warp in the outer part of its gaseous disk l lCarignan & Purtonl 19981 ). 
Nevertheless, the solid and dashed lines achieve excellent fits up 
to 6 kpc and the dashed line uses an acceptable distance and gas 
scale-height of only 0.65 kpc. The dotted line uses a razor thin gas 
disk a nd is a poor fit to t he inn er curve and poorer still to the outer 
curve. IWalter & Brinks! i ll 9991) claim that dwarf galaxies (such as 
DDO 154) have thicker H I disks than more massive spiral galaxies 
due to a lower restoring force from weaker gravity. In order to be 
a good fit in MOND, the outer warp should be responsible for the 
anomalously low rotation speeds between 7 and 8.5 kpc. 



5 RESULTS 

In Fig [5] we plot the minimum reduced \ 2 for each free pa- 
rameter as a function of that parameter. Specifically, this is the min- 
imum \ 2 achievable with that parameter fixed and every other pa- 
rameter free. The distances and M/L are limited according to the 
plotted range for each galaxy. 

Each galaxy in Fig|5]has a different line colour. In the top left 
hand panel, we plot reduced \ 2 f° r eacn galaxy against distance. 
The fitted distance is normalised to the measured distance and sim- 
ilarly the fitted M/L is normalised to the M/L predicted by the 
diet-Salpeter IMF. The la error on distance varies from galaxy to 
galaxy and can be found in Table Q] For M/L, the Kroupa IMF is 
typically 70% of the diet-Salpeter value. 

Clearly two galaxies (NGC 2903 and NGC 3521) prefer sub- 
stantially lower distances than the measurements suggest, whereas 
the other three galaxies have their minima closer to the measured 
distance. It is worth baring in mind that NGC 2903 and NGC 3521 
have the largest uncertainties on their distances, 25% and 30% re- 
spectively. 

The M/L plots (top right and middle right) have only four 
lines since we fix the M/L of DDO 154 to the diet-Salpeter value 
and do not vary it, nor do we vary its stellar scale-height. The con- 
cern here is that the curve for NGC 2903 rises sharply towards large 
X 2 when typical M/L are tried. 

The gas scale-height is only significant for two galaxies. Ob- 
viously DDO 154 is our only gas dominated galaxy and 95% 
of its mass is gaseous. It prefers as large a gas scale-height as 
possible and the \ r ' ses sharply below 0.6 kpc. On the other 
hand, NGC 2403 prefers a thin gas disk below 0.4 kpc in scale- 
height. The other three galaxies prefer gas scale-heights larger than 
0.5 kpc, but no significant gain is achieved for larger values than 
this. 

The stellar scale-height is strongly constrained by the rotation 
curves. All four stellar dominated galaxies show a clear minimum 
between 0.25 and 0.55 kpc and the quality of the fits would be 
significantly reduced if they were forced to be razor thin or 1 kpc 
in extent. 

In Figs 1 61101 we plot various fits to the five rotation curves 
as well as three contour plots (DDO 154 only has one) for the 
three combinations of the three free parameters (gas scale-height 
has been omitted). 



5.2 NGC 3621 

The % 2 minimum for NGC 3621 is obvious in the contour 
plots of Fig [7] and puts strong constraints on all three parameters. 
The fits, whose parameters are given in table Q] at the top left hand 
panel show the best fit (solid line) and another fit with slightly 
larger distance (dashed line), more in keeping with the measured 
distance. Both are excellent fits. The dotted line shows the best fit 
model when the best fit stellar scale-height is replaced with a razor 
thin stellar disk and is a very poor fit to the central rotation curve. 
The best fit parameters are all within acceptable bounds and this is 
an excellent fit for MOND. 

5.3 NGC 2903 

The contour plots in Fig [8] show that in order to have a satis- 
factory fit, regardless of distance and stellar scale-height, the M/L 
must be more than twice the predicted value. If this is the case, 
then an excellent fit can be achieved as shown in the top left panel. 
On the other hand, using a M/L of only 50% larger than the pre- 
dicted value gives a poor fit, as the dotted line attests. Therefore, 
NGC 2903 is a problem for MOND unless a reason can be found 
why it should have a M/L that is twice the diet-Salpeter prediction. 
If such a reason was found, it would be well fit by MOND. 

5.4 NGC 2403 

In Fig [9] the contour plot of distance against M/L shows the 
requirement for a good fit to have either an excessively large M/L 
or distance. This is aggravated by the tight error on the distance 
of a mere 8% and the relatively low predicted M/L for the given 
stellar population. For instance, the predicted 3.6 /im diet-Salpeter 
IMF M/L of NGC 2403 is 0.39, whereas for NGC 2903, 3521 and 
3621 it is 0.61, 0.73 and 0.59 respectively. 

In the top left hand panel of Fig [9] we fit four curves to the 
observed rotation curve. The solid and dashed curves are the best 
fit models, which respectively use a > 2a discrepant fitted distance 
and regular M/ L, and a more regular distance but more discrepant 
M/L. The other two lines are a model with the same parameters as 
the best fit, but using a razor thin stellar disk (dotted line), and one 
with a regular distance and M/L (dot-dashed line). Both of these 
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are unsatisfactory fits. As with NGC 2903, this galaxy NGC 2403 
can only be consistent with MOND if the cause of its higher than 
expected M / L is found, however, this is a much less significant 
problem than NGC 2903 because even with the required M/L for 
a good fit, it would still have a typical M/L compared to other 
galaxies. 

5.5 NGC 3521 

For NGC 3521, the preference is for a much lower distance 
than measured, which can be seen in the contour plot of Fig [10] 
This is not necessarily out of the question since the error on its 
distance is 30%, but the x 2 minimum is actually further from the 
mean than this. The parameters of the two models plotted against 
the observed rotation curve in the top left panel are given in tableQ] 
and correspond to a distance more than la lower than the mean and 
roughly 0.5a lower than the mean. Both require sensible M/L and 
the stellar scale-height is constrained to be close to 0.25 kpc. The 
fit with the less discrepant distance (dashed line) is a better fit to 
the inner 13 kpc and the discontinuity in measured rotation speed 
at this radius makes it difficult for any model to achieve a good fit 
to both the inner and outer curve. In general, this galaxy is well fit 
by MOND. 



5.6 Stellar scale-heights 

A somewhat surprising result is the tight constraint placed on 
the values of the stellar scale-height for the four stellar dominated 
galaxies and the gas scale-height of DDO 154. From the contour 
plots of Figs[6]to[T0]we see the relative lack of freedom in the stel- 
lar scale-heights and this gives us the opportunity to put a further 
constraint on the model by checking if the required scale-heights 
are consistent with the typical stellar vertical velocity dispersions 
of other spiral galaxies. 

Making the assumption that the stellar velocity dispersions are 
isothermal with height above the disk we can use the following 
equation to solve for the stellar vertical velocity dispersion at any 
radius 



d$(R,z) p*(R,z) 



dz 



d zP *(R,z) 



(5) 



A similar test could be used to constrain the gas scale-heights, but 
we don't follow that route. The right hand side is comprised fully 
of quantities that are known in the rotation curve fitting code. As 
can be seen in Fi d 1 1 1 the vertical velocity dispersions of all the 
galaxies have a similar shape. Generally, for the best fit models 
(solid and dashed lines) the trend is to have a large central vertical 
velocity dispersion a z (0) of between 20 and 35kms _1 and this 
drops to between 4 and 8 km s~ 1 by 15 kpc. The dotted lines, which 
are for very thin disks of only 0.05 kpc scale-height, have much 
lower velocity dispersions - with NGC 2403 and NGC 3621 being 
cent rally just 8 and 1 kms -1 respectively. 

iBottemal d 19931) studied the stellar velocity dispersions of a 
sample of 12 spiral galaxies of varying blue band magnitude in the 
range Mb =-18.76 to -22.22. They found a trend for the vertical 
stellar velocity dispersions to scale with decreasing galaxy magni- 
tude (increasing luminosity). The three galaxies with Mb > —20, 
for which only NGC 3521 from our sample does not apply, were 
NGC 3 198, 3938 and 6503 and they had central stellar velocity dis- 
persions of 45 ± 5 kms - , 32 ± 13 kms -1 and 30 ± 7 kms" 1 
respectively, although NGC 3198 and NGC 3938 are inclined. The 



stellar velocity dispersions dropped to between 10 and 20 kms -1 
after 2 or 3 scale-heights, which is similar to our models. The larger 
galaxies with Mb < —21, which applies to NGC 3521, had cen- 
tral stellar velocity dispersions between 50 and 120 kms -1 . This 
is larger than the central stellar velocity dispersion of NGC 3521. 
Even increasing the stellar scale-height to 1 kpc can only increase 
the central vertical stellar velocity dispersion to 6 kms -1 and thi s 
remains on the low end of the sample studied by |Bottemj ( ll993h . 
Perhaps the increased velocity dispersions with luminosity are also 
linked to the increased prevalence of bulges with luminosity and 
NGC 352 l's lack of a significant spheroid bucks this trend. 

The stellar and gas scale-heights we use in the fits are 
given in Table [J alo ng with the stellar scale-heights suggested by 
Ide Blok et al.l d2008l ). The suggested stellar scale-heights are sim- 
ply one-fifth of the scale-lengths and for NGC 3521 and NGC 3621 
respectively are 1.9 and 1.5 kpc. These values are several times 
larger than the fitted scale-heights of 0.25 and 0.5 kpc respec- 
tively. If these suggested scale-heights were enforced, it would not 
be possible to achieve a good fit to the rotation curves. There- 
fore, it is vital to have a separate measurement of the stellar ve- 
locity disper sion to put an orthog onal limit on the scale-heights, 
as shown bvlPuglielli et al.l fcOld) in t he cas e of NGC 6503 and 
iBershadv et alj d201ll) . lBershadv etal] J201ll) used the Disk Mass 
Survey to demonstrate that the magnitude of vertical velocity dis- 
persions led to the inference that galaxy disks must be submaximal. 
Maximum disk simply means the fit to the rotation curve with the 
highest M/L possible and the least amount of dark matter, such 
that all the rotation velocity is attributable to the disk in the central 
part. It would be interesting to see if this holds in MOND for which 
the disk, by definition, mu st be maximal. 

iBaneriee et al. I J201 J) made an analysis of the scale-heights of 
the gas dominated THINGS galaxies in Newtonian gravity and it 
is essential that a similar review is made of the same galaxies in 
MOND. 



5.7 The external field effect 

External gravitational fields (see forthcoming paper), usually 
the result of nearby galaxies or clusters, can cause a suppression of 
the boost to gravity due to MOND. This is particularly likely for 
low surface brightness galaxies. An external gravitational field has 
a similar effect to decreasing the M/L and increasing the scale- 
height, but they are not degenerate given the gaseous mass content 
and the radial dependence of the circular velocity on scale-height. 
It is an oft-quoted solution to unknown problems in MOND, but it 
is important to emphasise that it would have no beneficial influence 
on NGC 2403 or NGC 2903 since it would simply impose a larger 
M/L, making the situation worse. 



6 CONCLUSION 

Here we have introduced an N-body co de that solves th e mod- 
ified Poisson equation of QUMOND (see lMilgronj[201fj|) . With 
it, we fitted the rotation curves of five spiral galaxies from the 
THINGS survey dWalter et alj|2008l) , using N-body realisations of 
the stars and the gas in each galaxy fixed by their surface densities. 
We allowed the distance, mass-to-light ratio and both the scale- 
heights of the stellar and gaseous disks to be free parameters, with 
priors set depending on observational constraints. 

We discovered that our best fits were excellent matches to the 
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Galaxy 


Measured distance 


Fitted distance 


Stellar (Gas) Mass 


Predicted M/L 


Fitted M/L 


Suggested z* 


Fitted z t (z g ) 


Line-type 


Measured distance 


dS M/L 




Mpc 




10 9 M Q 


3.6^m d-S (Kr) 


kpc 


kpc 




DDO 154-a 


4.30 ± 1.07 


1.085 


0.026 (0.468) 


0.32 (0.23) 


1.0 


0.2 


0.2(1.5) 


Solid 


DDO 154-b 




0.96 






1.0 




0.2 (0.65) 


Dashed 


DDO 154-c 




0.9 






1.0 




0.2 (0.05) 


Dotted 


NGC 3621 -a 


6.64 ± 0.70 


0.91 


19.3 (9.58) 


0.59 (0.42) 


0.88 


1.5 


0.5 (1.0) 


Solid 


NGC 362 1-b 




0.97 






0.76 




0.5 (1.0) 


Dashed 


NGC 362 1-c 




0.91 






0.88 




0.05 (1.0) 


Dotted 


NGC 2903-a 


8.90 ± 2.20 


0.8 


16.2(6.6) 


0.61 (0.43) 


4.6 


0.5 


0.54 (1.0) 


Solid 


NGC 2903-b 




1.1 






2.5 




0.3 (1.0) 


Dashed 


NGC 2903-c 




1.15 






1.5 




0.2(1.0) 


Dotted 


NGC 2403-a 


3.47 ± 0.29 


1.19 


5.13 (3.82) 


0.39 (0.26) 


1.15 


0.4 


0.4 (0.1) 


Solid 


NGC 2403-b 




1.09 






1.45 




0.5 (0.1) 


Dashed 


NGC 2403-c 




1.19 






1.15 




0.05 (0.1) 


Dotted 


NGC 2403-d 




1.05 






1.2 




0.6(0.1) 


Dot-Dashed 


NGC 3521-a 


10.7 ± 3.20 


0.63 


125.5 (13.0) 


0.73 (0.52) 


1.1 


1.9 


0.25 (1.0) 


Solid 


NGC 3521-b 




0.82 






0.8 




0.25 (1.0) 


Dashed 



Table 1. Here we show the various parameters corresponding to our fitted models that are plotted in Figs [6lto ll0t . For each galaxy we give the measured 
parameters from the literature and our fitted values. The measured distances come from Walter et al. (2008), who took their value from Freedman et al. (2001), 
except NGC 2403 which comes from the more recent study of Vinko et al. (2006). The masses for the stellar and gaseous components come from de Blok et al. 
(2008) and are correct at the mean of the measured distance and with use of the diet-Salpeter initial mass function. The predicted mass-to-light ratios (M/L) 
also come from de Blok et al. (2008), for which "d-S" and "Kr" correspond to diet-Salpeter and Kroupa initial mass functions respectively. The suggested 
stellar scale-heights are merely one-fifth of the radial scale-height. There is no suggested scale-height for the gaseous distribution. For each galaxy, every 
separate fit is given a different a letter - a, b, c or d - and these fits correspond to different linetypes in their corresponding figures. Our fitted distances and 
values for the M/L are normalised to the mean measured distance and the diet-Salpeter IMF respectively. The stellar scale-heights are well constrained, but 
the gas scale-heights less so. For DDO 154 we fixed the stellar parameters since they had little influence. 



rotation curves of all five galaxies, except for some minor discrep- 
ancies at locations of publicised uncertainty in the observations. We 
displayed contour plots of reduced y 2 for the various free param- 
eters of all galaxies and found that two galaxies, NGC 2403 and 
NGC 2903, could only make satisfactory fits if their mass-to-light 
ratios were larger than the predicted values. The required increase 
in M/L for NGC 2403 is moderate, but for NGC 2903 it is con- 
siderabl e. Interestingly , the da rk matter fits to these rotation curves 
given bv ldeBloketaU ( l2008h . which can be seen in their Tables 5 
and 6, also require larger M/L. For NGC 2403, a value 50% larger 
than the predicted diet-Salpeter value and for NGC 2903 more than 
double was required, which is exactly what we have found. 

If this were an isolated incident, it could have less relevance, 
but there is a growing preference for stellar dominated galaxies 
to require much larger M/Ls than predicted by stellar population 
models from typical initial stellar mass functions. Large M/Ls are 
also required for certain dwarf spheroidal galaxies surround ing the 
Milky Way, in particu lar Carina, Sextans and Draco (see I Angus! 
l200alSerra et al.l20ld) and a large frac t ion of the early type galax- 
ies studied bv lSanders & Noordermeeil d2007f) . The MOND fits to 
the dynamics are in general still excellent, but a solution must be 
found to explain why certain galaxies can have M/L described by 
a Kroupa or diet-Salpeter initial mass function and other galaxies 
need one up to twice as large. Both the MOND and dark matter 
fits to the rotation curves of NGC 2403 and NGC 2903 suggest 
there is a problem with the stellar population synthesis models of 
these galaxies. Furthermore, the distance required for NGC 3521 is 
considerably lower than the mean, but this is acceptable due to the 
large uncertainty. A revised distance with tighted error bars could 
be very revealing. These issues need to be resolved before we can 
say that MOND provides good fits to all five galaxies with reason- 
able parameters. 

The surprising result was that the MOND fits put a tight con- 
straint on the stellar scale-heights. The best fit scale-heights of all 
four stellar dominated galaxies were found to be between 0.25 kpc 



and 0.55 kpc and were strongly constrained to be larger than and 
less than 1 kpc. For two galaxies, NGC 2403 and 2903, the best 
fit stellar scale-heights were very close to the typical scale-height 
used, which is one fifth of the radial scale-length. On the other 
hand, NGC 3521 and NGC 3621 have best fit stellar scale-heights 
that are far lower than the suggested scale-height insofar as if the 
scale-height was fixed to this suggested value, a good fit to the inner 
data points would not be achievable. Similarly, the only gas domi- 
nated galaxy we studied, DDO 154, can achieve a good fit only if 
its gas scale-height is larger than ~0.6 kpc. In fact the quality of 
the fit improves with increasing scale-height. 

Clearly, the use of both stellar and gas scale-heights as free 
parameters in MOND fits to galaxy rotation curves is imperative. 
Line of sight velocity dispersions of both components in a sam- 
ple of similar face on galaxies should be used as a sanity check on 
their values since in MOND the rotation curves and vertical veloc- 
ity dispersions of both components are bound together to the mass 
distribution. 
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the modified Poisson equation of QUMOND (solid black, Eq[2} and the 
algebraic relation (solid red, Eq[5J for a razor thin exponential disk. The 
dashed curve shows the rotation speed from the same exponential disk, us- 
ing our code, but with a vertical scale-height of 1 kpc that is described by a 
sech 2 distribution. 
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continuities are caused by moving from one grid resolution to another. 
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Figure 5. Shows the curves of minimum reduced \ 2 as a function of the fitted parameters of each galaxy. In the top left panel is the fitted distance of each 
galaxy normalised to the observed distance, and the top right is the fitted M/ L normalised to the M/L using the diet Salpeter IMF (the figure below this is the 
same as above but with a larger plotted range). We also plot the two scale-heights. A different coloured line is used for each galaxy and a legend emphasises 
this in each panel. 
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Figure 6. In the left hand panel we plot the rotation curve of DDO 154 as measured by de Blok et al. (2008) along with our MOND fits. The curves are 
described in table[T] giving their line-types and model parameters. The right hand panel shows contours of reduced \ 2 for the two free parameters: the distance 
normalised by the mean measured distance and the gas scale-height. The contour levels are defined in the panel and the vertical dashed line defines the la 
error on the measured distance. 
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Figure 4. Shows the surface densities of the stellar (solid) and gaseous 
(dashed) components for NGC 3621. From de Blok et al. (2008). 
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Figure 7. In the top left hand panel we plot the rotation curve of NGC 3621 as measured by de Blok et al. (2008) along with our MOND fits. The curves are 
described in table[T] giving their line-types and model parameters. The other three panels show contours of reduced \ 2 f° r the three combinations of three free 
parameters, these are the fitted distance and M / L normalised to the measured distance and diet-Salpeter IMF respectively and the fitted stellar scale height. 
The contour levels are defined in the top right panel. The vertical dashed line in the contour plots with distance defines the la error on the measured distance. 
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Figure 8. As per Fig|7] but for NGC 2903. 
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Figure 9. As per Fig|7] but for NGC 2403. 
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Figure 10. As per Fig|7] but for NGC 3521. 
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Figure 11. Stellar vertical velocity dispersions as functions of radius for the four stellar dominated galaxies. The linetypes for each galaxy are defined such 
that the solid, dotted and dashed lines always corresponds to models a, b and c respectively as set out in table[T] 



